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ABSTRACT 


Linear  search  algorithms  are  developed  for  use  when  minimizing 
logarithmic  barrier  functions,  whose  one-dimensional  taLwior  is  in 
general  modeled  poorly  by  the  low-order  polynomial  approximations 
of  standard  linear  search  procedures.  The  new  methods  are  based  on 
special  approximating  functions  with  a logarithmic  singularity,  and 
are  designed  to  utilize  the  same  information  as  procedures  based  on 
quadratic  or  cubic  polynomials.  Although  the  parameters  of  the  special 
approximating  functions  depend  nonlinearly  on  the  available  data,  the 
determination  of  the  parameters  requires  little  additional  woric  in 
comparison  with  polynomial  fits.  Use  of  the  special  approximating 
functions  has  led  to  a significant  improvement  in  efficiency  when 
minimizing  logarithmic  barrier  functions,  where  efficiency  is  measured 
by  the  number  of  function  (or  function  and  gradient)  evaluations 
required  for  termination  of  each  linear  search. 


ACCESSION  for 


1/ 


1.  Introduction 


An  essential  part  of  many  algorithms  for  minimizing  a function 
of  several  variables  is  the  determination  of  a positive  step  to  be 
taken  from  an  initial  point  along  a given  search  direction.  At  the 
new  point,  the  values  of  the  function  and  possibly  its  derivatives 
are  required  to  satisfy  some  specified  conditions  with  respect  to 
these  quantities  at  the  original  point.  For  example,  in  theoretical 
presentations  of  such  algorithms,  it  is  often  stated  that  the  step 
taken  should  minimize  the  function  along  the  direction  of  search) 
the  desired  step  in  this  case  will  be  denoted  by  Of*,  and  is  defined 
by: 

F(x  + o*p)  = min  F(x  + ap), 

a 

where  x is  the  initial  point,  p is  the  search  direction,  and  F(x) 
is  the  function  to  be  minimized. 

However,  in  most  applications  the  conditions  that  define  an 
acceptable  step  are  less  stringent,  and  will  be  satisfied  by  an  approx- 
imation to  CP  (see  Osborne,  1972)  Gill  and  Murray,  197*0*  An  efficient 
method  for  obtaining  such  a step  can  be  based  on-  an  iterative  procedure 
to  determine  Of*,  by  generating  iterates  only  until  the  required  con- 
ditions hold.  Thus,  whatever  the  criteria  to  be  satisfied,  the 
underlying  iterative  process  remains  the  same)  what  varies  is  the 
number  of  elements  in  the  sequence  that  are  actually  computed. 


The  discussion  to  follow  will  be  concerned  only  with  the  procedure 
by  which  successive  estimates  of  the  desired  step  length  are  generated, 
and  not  with  the  termination  criteria.  The  terms  "linear  search" 
and  "one-dimensional  minimization"  will  be  used  to  denote  the  iterative 
procedure  common  to  all  such  algorithms,  and  do  not  imply  that  a 
close  approximation  to  Of*  is  to  be  obtained. 

The  linear  search  procedure  is  of  crucial  importance  to  the 
success  of  most  computational  algorithms  for  unconstrained  minimization) 
consequently,  there  has  been  considerable  work  on  devising  efficient 
and  reliable  methods  for  one -dimensional  minimization,  and  on  the 
careful  implementation  of  such  methods  as  computer  routines  (see  Brent, 
1973)  Gill  and  Murray,  197^) • These  general-purpose  linear  search 
algorithms  combine  successive  low-order  polynomial  interpolation  or 
extrapolation  with  safeguards  to  ensure  convergence  and  numerical 
stability.  The  performance  of  a linear  search  procedure  is  normally 
measured  in  terms  of  the  number  of  function,  or  function  and  gradienx. 
evaluations  required  to  locate  a point  that  satisfies  the  prescribed 
termination  criteria.  When  executing  a linear  search  with  respect 
to  a particular  nonlinear  function,  the  efficiency  of  the  algorithm 
will  thus  depend  on  the  adequacy  of  the  polynomial  approximation  to 
the  given  function.  For  minimizing  general  functions  with  no  a priori 
information  about  the  nonlinearities,  linear  search  procedures  based 
on  polynomial  approximation  may  be  the  most  efficient,  although 
inevitably  the  non-polynomial-like  behavior  of  some  functions  will 
necessitate  a large  number  of  successive  polynomial  fits  (and,  hence, 
evaluations  of  the  function). 

? 


This  paper  will  be  concerned  with  the  design  of  linear  search 
algorithms  for  a particular  class  of  nonlinear  functions — logarithmic 
barrier  functions.  A linear  search  with  respect  to  these  functions 
is  required  in  several  algorithms  for  constrained  minimization,  and 
yet  their  behavior  is  in  general  modeled  poorly  by  low-oivjr  polynomials. 


2.  The  Logarithmic  Barrier  Pi 


Barrier  function  methods  were  first  introduced  as  a technique 


for  transforming  a nonlinearly  constrained  minimization  problem  into 


a sequence  of  unconstrained  problems,  where  feasibility  with  respect 
to  the  problem  constraints  is  maintained  throughout  (a  full  description 
is  given  in  Fiocco  and  McCormick,  1968) . The  only  barrier  function 
to  be  considered  here  is  the  logarithmic  barrier  function,  defined 


as  follows:  corresponding  to  the  nonlinearly  constrained  minimization 


problem 


minimize  F(x) , x € En 

(2.1) 

subject  to  c.(x)  > 0,  i - 1,  2,  •••,  £, 

where  F(x)  and  (c^(x)}  are  prescribed  nonlinear  functions,  the 
logarithmic  barrier  function  is  given  by: 

l 

(2.2)  B(x,  r)  =?  F(x)  - r £ Hn(c.(x)), 

i -1 

where  the  positive  scalar  r is  termed  the  "barrier  parameter." 


Interest  in  barrier  function  methods  has  waned  in  recent- 


years,  due  to  the  development  of  more  promising  algorithms,  many 
of  which  are  based  on  extensions  to  the  quadratic  penalty  function 
(see  Fletcher,  197^j  for  a discussion  of  penalty  functions,  see  Fiacco 
and  McCormick,  1968) . However,  nearly  all  the  more  recent  algorithms 
are  non-feasible  point  methods,  and  hence  are  unsatisfactory  for 
problems  in  which  the  objective  and/or  constraint  functions  are  ill- 
defined  or  undefined  outside  the  feasible  region.  In  addition,  non- 
feasible  point  methods  may  have  other  disadvantages.  For  example, 
the  transformations  employed  may  introduce  spurious  non-feasible 
solutions j furthermore,  if  such  an  algorithm  is  terminated  prematurely, 
it  often  fails  to  provide  a useable  solution. 

Many  of  the  algorithmic  developments  based  on  the  quadratic 
penalty  function  can  be  mirrored  by  utilizing  the  logarithmic  barrier 
function.  Two  such  suggestions  have  been  made  by  Osborne  (1972)  and 
Murray  and  Wright  (1976).  For  feasible  point  algorithms,  the  logari- 
thmic barrier  function,  which  has  the  virtue  of  maintaining  feasibility, 
can  often  serve  as  a convenient  "merit  function"  r0r  a linear  search, 
regardless  of  how  the  search  direction  is  obtained.  Therefore,  it 
is  worthwhile  to  consider  how  to  exploit  the  special  properties  of 
the  logarithmic  barrier  function  in  the  design  of  a linear  search 
algorithm. 

The  inadequacy  of  linear  search  procedures  based  on  polynomial 
fits  for  minimizing  barrier  functions  has  been  discussed  by  Fletcher 
and  McCann  (1969),  Murray  (1969),  Lasdon  et  al.  (1975),  and  Ryan  (197*0  • 

h 


The  defining  characteristic  of  a barrier  f'i:vt.i:jn  i.-  uinguiarity 
at  the  boundary  of  the  feasible  region;  hi  nee  polynomials  provide 
a poor1  approximation  r.o  a function  with  a singularity,  the  general 
approach  has  been  to  suggest  alternative  approximating  functions 
with  the  same  kind  of  singularity  ns  the  barrier  function  The 
minimum  of  the  special  appro” imating  function  can  be  taken  as  an 
estimate  of  the  minimum  of  the  barrier  function,  to  be  used  within 
the  linear  search  procedure  In  exactly  the  same  way  as  the  estimated 
minima  of  polynomial  approximations.  However,  it  may  not  be  possible 
to  determine  explicitly  the  coefficients  and/or  the  minimum  of  the 
new  approximating  function,  in  contrast  to  these  calculations  for 
polynomials. 

5.  Fitting  of  Special  Functions 
5.1.  Dlscuasion 

The  special  functions  w ot  considered  are  designed  to  contain 
a single  logarithmic  singularity.  This  restriction  to  a particular 
barrier  function  contrasts  with  the  general  approach  taken  by  Lasdon 
et  si.  (197?),  where  a special  function  is  developed  by  the  applica- 
tion of  the  form  of  the  barrier  function  to  linearized  approximations 
to  the  objective  function  and  the  constraints.  The  approach  taken 
here  allows  a much  simpler  solution  to  the  problem,  where  the  known 
form  of  the  singularity  is  directly  exploited.  A single  function  is 
used  to  approximate  the  behavior  of  the  barrier  function  along  the 
search  direction,  and  it  is  not  necessary  to  make  separate  approxima- 
tions to  the  objective  function  and  the  constraint  functions. 
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The  logarithmic  barrier  function  can  only  at 

feasible  points,  and  Is  undefined  beyond  *.he  fir.-*,  oinrilarity 
along  a particular  direction.  Therefore,  only  feasible  points  are 
considered  in  fitting  the  approximating  fune’iotis,  since  the  lata 
to  be  used  will  be  values  of  the  barrier  function  and  -to  gradient. 

It  will  be  assumed  that  we  seek  an  approximation  to  the 
minimum  along  the  search  direction  of 

£ 

B(x,  r)  * F - r E in(c .) , 
i-1 

where  the  barrier  parameter,  r.  Is  known  and  fixed  throughout  the 
linear  search.  The  special  functions  to  be  considered  are  the  follow- 
ing, where  the  parameter  8 represents  the  one-dimensional  variation 
along  the  search  direction: 

(a)  8 linear  function  plus  a logarithmic  singularity,  of  the  form: 

fL(9)  = S + $6  - r /n(d  - 8)> 

(b)  a quadratic  function  plus  a logarithmic  singularity,  of  the  form: 

fn(8)  = a + b9  + c02  - r £n(<f  - 9). 

Thses  formulations  attribute  all  the  singular  behavior  to 
one  location  (8  = in  fact,  the  singularity  displayed  by  a 
barrier  function  depends  in  general  upon  the  first  constraint  zero 
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encountered  along  a gi\en  direction,  and  ieros  that  may  occur  beyond 
the  first  are  irrelevant.  This  behavior  becomes  especially  marked  for 
small  values  of  the  barrier  parameter,  when  the  constraint  functions 
bounded  away  from  zero  have  almost  no  influence  on  the  barrier 
function,  and  the  nearest  constraint  has  influence  only  very  close 
to  its  zero  along  the  given  direction.  A smooth  function  with  a 
suitable  damped  singularity  should,  therefore,  be  an  excellent  model 
of  a barrier  function  close  to  the  boundary  of  the  feasible  region. 


Plus  Logarithmic  fiingularll 


The  special  function  to  be  fitted  is  of  the  form: 


fT  = a + - r in(cf  - 9). 


Differentiating  f^  with  respect  to  9 gives: 


3 - 9 


'L  " (5  - 9)2 


Since  the  barrier  function  is  decreasing  at  0=0,  this  same  condi- 
tion will  be  required  of  f^,  so  that 


f-(°) 


b + ? < 0. 

Q 


The  expression  for  is  equal  to  the  value  of  f J (C)  times 

(d/£) , and  since  (£/b)  and  f£(0)  ore  negative,  9*  oust  he 
positive.  f£  is  everywhere  positive,  so  that  the  stationary'  point 
must  be  a minimum  of  f^. 

In  order  to  specify  the  function  f j , three  independent  pieces 
of  information  are  required  about  the  behavior  of  the  function  to 
be  fitted  along  the  search  direction;  f^  will  therefore  be  used  in 
circumstances  where  a parabolic  fit  would  normally  be  carried  out  to 
minimize  a general  function.  Figure  1 illustrates  a barrier  function 
and  its  approximations  by  the  function  f^  and  a parabola,  l>oth  using 
the  same  data.  It  is  clear  that  the  special  function  gives  a more 
accurate  prediction  than  the  parabola  of  the  minimum  of  the  barrier 
function. 
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r = .005 


Figure  1 
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The  special  function  fT  will  be  fitted  *.»ith  the  same  sets 
of  data  used  by  a typical  parabolic  line  search,  namely:  (1)  one 

function  value,  the  corresponding  projected  gradient,  and  a second 
function  value;  (2)  three  function  values.  In  the  following  discus- 
sion, f,  will  denote  the  function  value  at  9.,  and  g.  is  the 

i i i 

corresponding  value  of  the  projected  gradient. 

Case  JL  -Two  function  values,  one  gradient 

The  three  unknown  parameters  of  f^  — a,  b,  and 
solved  for  in  terms  of  the  known  values.  Assuming  that 
the  equations  specifying  f^  are: 

(A)  f^  * a + b0^  - r in(d  - 0^) 

(b)  ex  ■ * + 

(c)  r2  = a + b02  - r £n(cf  - 0g) . 

There  is  no  loss  of  generality  in  assuming  that  0^  = 0,  so 

that  from  (A)  and  (B)  we  obtain: 

a = + r ln(d; 


d — must  oe 


b2  > V 


and 


b 


10 


I The  problem  to  be  solved  then  becomes  that  of*  finding  a solution  d 

? to  satisfy 

i 


Because  ©£  is  feasible,  the  value  of  the  singularity,  £, 
must  lie  in  the  interval  (0£,  *).  The  function  •:  ^(d)  has  the 

following  properties: 


-*  - » as  d -»  + * 


4^  -*  as  d -*  • * 


> 0 


for  all  finite  dj 


(Jj^  -♦  so 


as  d -+  + } 


-*  0 


as  d -*  » ) 


< 0 for  all  finite  d. 


These  relations  imply  that  there  is  a unique  zero  of  $>  in  *) » 

provided  that  k1  > 0,  l.e.,  (f£  - fj)/^  > gi*  ’ntis  r*Juirem*mt 
means  that  the  function  value  at  0£  must  be  larger  than  the  linear 
approximation  at  9^  would  have  predicted,  i.e.,  f£  must  lie  in 
the  shaded  region  of  Figure  2. 


Although  the  numerical  solution  of  the  equation  4>^(d)  » 0 

A 

would  yield  the  desired  value  4,  t is  an  ill-behaved  function, 
unsuitable  for  the  usual  zero-finding  techniques  such  as  Newton's 


method.  Figure  5 illustrates  the  behavior  of  <1^. 


An  attempt  to  apply  Newton’s  method  in  Region  I would  cause 


very  slow  convergence  to  the  zero,  although  all  estimates  would 
undershoot  and  hence  could  not  diverge.  In  Region  II,  Njwton's 
method  might  easily  yield  an  estimate  to  the  left  of  Sg,  and  safe- 
guards would  need  to  be  incorporated  to  prevent  divergence. 

To  avoid  these  difficulties,  the  problem  of  solving  ^(d)  = 0 
c&n  be  transformed  into  an  equivalent  problem  that  is  easy  to  solve. 
If  we  define  a new  variable  y * Sg/d , so  that  0 < y < 1 for  the 
admissible  range  of  d,  the  equation  ®^(d)  ^ 0 can  be  written  in 
terms  of  y as 

ln(l  - y)  - -y  - kx  . 

Taking  exppnentials  of  both  sides  gives  the  equation 


and  the  value  ^ which  satisfies  this  equation  is  a zero  of  a new 
function. 


♦x  (y) 


Differentiate  with  respect  to  y gives: 


and  thus  and  V"  are  negative  for  all  y In  the  interval 

(0,  l) , if  > 0.  Furthermore, 

“kl 

»1(0)  - 1 - e >0, 

-1  “kl 

f1(l)  - -e  e < 0, 

so  that  has  a unique  aero  y in  (0,1)  if  > 0.  Figures  t 
and  5 illustrate  the  behavior  of  for  two  values  of  k^.  Even 
for  k^  small,  when  the  nonlinear  portion  of  is  more  significant, 

is  very  close  to  linear  in  (0,  l). 

Because  is  well-behaved,  Newton’s  method  will  converge 

to  y very  rapidly.  However,  even  further  advantage  can  be  taken 
of  the  form  of  f ^ by  noting  that  y satisfies  the  relation 

(3.2)  kx  = -y  - £n(l  - y) . 

Since  the  right-hand  side  of  (?.2)  depends  only  on  the  parameter 
y,  and  is  independent  of  the  problem  data,  the  vslues  of  the  function 
cn  the  right-hand  side  can  be  tabulated  for  a set  of  y in  (0,  l). 

For  a particular  value  of  k^,  table  lookup  could  then  be  used  to 
determine  a highly  accurate  estimate  of  y.  With  a sufficiently 
large  table,  no  iteration  would  be  necessary  to  locate  the  solution 
to  any  desired  accuracy^  but  since  Newt on* s method  converges  rapidly, 


1C 


^,(y),  k,  =1.4 
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a small  table  plus  a single  iteration  will  locate  y (and  d)  to 

the  required  accuracy  (e.g.,  seven  values  are  sufficient  for  an 

-5 

accuracy  of  10  ) . 

Case  2.  Three  Function  Values 

Assuming  that  9^  > 9^  > 6^,  the  three  equations  specifying 


the  parameter:; 

•*'»  h, 

arid 

d of  f^  are: 

(A) 

f 

Ll 

A 

- a 

+ bO  - r in(d  - 0 ) 

<B) 

*2 

A 

■t  a 

+ £<?,.,  - r in(5  - 0 ) 

r.'  c 

(C) 

f5 

- a 

+ b9  - r £n(d  - 9 ) 
2 ✓ 

There  Is  no  lose  of  generality  in  assuming  that  0^  = 0. 


Us! rig  equfctlu' 


* .*  el itai note  o,  we  obtain: 


(A  nul  < u) 


These  two  equations  can  then  be  manipulated  to  eliminate  b,  yielding: 


(5.3) 


k 444  k -M ♦ i ( 


- f,  - f. 


a relationship  which  must  be  satisfied  by  &>  the  location  of  the  . 
singularity. 

Consider  (3*3)  as  a nonlinear  equation  in  terms  of  the  function 


•2(d) 


, /a  - 0 


‘)  + V 


where 


"i(: 


f,  - f,  f„  - f. 


The  problem  to  be  solved  then  becomes  that  of  finding  a solution  d 
to  satisfy  <&2(d)  = 0. 

Because  ©g  and  9 ^ are  feasible,  the  value  of  d must  lie 
in  the  interval  (0  , ») . The  function  $g  has  the  following  proper- 
ties : 

®g  -*  - “ as  d -♦  9^  + ) 

$2  -»  as  d -+  * ) 

$2  > 0 for  all  finite  dj 

•A  “* 00  as  d -»  0 + t 

3 

0 as  d -*  * ; 


r fo  if  ! 

L iiL 1 


*2  < 0 for  all  finite  d. 
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These  relations  imply  that  there  is  a unique  zero  of  <t>2 

in  (6y  80 ) if  kg  > 0.  This  requirement  has  a similar  interpretation 

to  Case  1,  i.e.,  the  configuration  of  f.,  f , and  f,  must  be  as 

x 2 5 

shown  in  Figure  6 , so  that  the  line  joining  f and  f1  must  lie 
above  the  line  joining  f£  and  f^ 

Because  of  its  highly  nonlinear  behavior,  4>2  is  unsuited  for 
application  of  the  usual  zero-finding  techniques.  The  behavior  of 
<X>2  is  illustrated  in  Figure  7. 

A 

To  transform  the  problem  of  finding  d to  solve  ^(d)  =»  0 
into  a computationally  manageable  form,  we  introduce  the  variables 
3 a d/0g,  so  that  3 > 1,  and  9^  = with  0^  > 1;  in  essence,  we 

consider  0g  to  be  unity  and  scale  0^  and  d accordingly.  If  we 
multiply  the  equation  $>2(d)  =0  by  and  substitute  the  variables 
d and  0,,  we  obtain: 

(5.10  to(i  -t)‘S  - l)  - V3  • 

let  the  variable  y = 0^/3,  so  that  0 < y < 1 for  all  admis- 
sible d.  Taking  exponentials  of  both  sides  of  (3.^)  gives: 


pm;Li  -ljpju  u"u,..  i-  m •, r . 3^7 . rr 


tprrM  WatV'-i 


• ■-  ■Tf*s»’’  .v-n- i-'ra* **•  •■^5:* 


where 


7 a > 1,  and  0 = 


-k„6 


2 3 


(so  that  0 < p < 1 for  k0  > 0) . 


Differentiating  with  respect  to  y,  we  obtain 


* -i  + e(i  - *)r  1 j 


7-2 


has  the  following  properties: 


*2(°)  = 1 - 0 > 0 


(y  = 0 corresponds  to  d-*  »)> 


«2d)  = -e(i  -y)r<° 


(y  = 1 corresponds  to  d -*  0^  +)  j 


*|(0)  - -1  + 0 < 0 


*|U)  - -1  + ?(l  - 


7-1 


< 0, 


and 


fi- 

fe 

If' 


w*r. 

> 

Hi' 


*"  < 0 


for  y in  the  interval  (0,  1). 
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These  conditions  imply  that  V n has  a unique  zero  in  (0,  l), 

which  could  be  located  by  Newton's  method.  H ^ is  a well-behaved 

function,  in  contrast  to  $0.  The  nonlinearity  in  ¥ results  from 

the  expression  (1  - y /y)^ , whore  y > 1,  and  is  the  quotient 

from  the  original  protkm.  The  value  of  y has  been  monitored  during 

several  runs  with  particular  barrier  functions;  it  never  exceeded 

2.0,  and  was  usually  in  the  range  (1.01,  1.5)  • Even  if  y were 

large,  the  expression  0 - y /y)^  approaches  e~^  as  y approaches 

",  and  hence  the  function  takes  on  the  form  of  ¥ in  Case  1. 

" 4 8 1 

Purthemore,  the  value  of  p(  = e J)  is  usually  small,  so  that 
the  nonlinear  component  of  is  even  less  significant. 


Figures  8 and  9 illustrate  the  behavior  of  f for  two  sets 
of  parameters  ((3,  y) . Note  that  even  for  large  3 and  y,  the 
nonlinearities  do  not  have  a significant  effect. 

It  is  nut,  computationally  convenient  to  obtain  a tabulation 
of  in  order  to  find  an  initial  estimate  of  y,  since  the  influences 

of  y and  y are  not  separable,  However,  Newton's  method  will 
converge  rapidly  from  a reasonable  starting  point  (for  example,  the 
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Differentiating  f^  with  respect  to  0 gives: 


?e  + 


(d  - ©y 


and 


2r 

.3  - e)3 


Since  the  barrier  function  is  decreasing  at  © =»  0,  this 
same  condition  will  be  required  of  f^,  so  that  f^(0)  < 0.  Hence, 
there  must  exist  a minimum  of  f^  in  the  interval  (0,  d),  because 
f_  -*  » as  © -*  5, 

w. 

A stationary  point  of  fn  will  occur  at  ©*  where 
f • (©*)  - 0,  so  that  0*  satisfies 

£ + 2°9*  + 3 e#  = o. 

This  relation  leads  to  a quadratic  equation  satisfied  by  0*, 

A 2 , A AA.  AA 

2c©  + (b  - 2cd)  0 - r - bd  = 0, 
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with  two  solutions  which  can  be  written  in  the  alternative  forms: 


(3.5)  e*  •- 


2c  d 


- b ± J(Y)  - 2cd ) -t  8c  ( r + b3)  , 
he 


or 


(3.6)  e*  = 


Pc<j  - b 


- 

4c 


fore  • 


We  now  examine  the  properties  of  these  two  solutions  in  order 
to  determine  which  will  qualify  as  a suitable  minimum,  i.e.,  lie  in 
the  interval  (0,  £) . Since  Cq(0)  < 0 and  d > 0,  the  expression 
£di  + r is  negative.  The  it*  fore,  the  sign  of  the  term  8c  (r  ■*  ££) 
is  opposite  to  the  sign  of  c.  Consider  two  cases: 


Case  A:  c > 0 (the  Quadratic  term  in  f has  a positive  coefficient). 

r\ 

Hie  quantity  (£>  + 2cd)<:  + 8rc  involved  in  fomula  (3.6)  for 
9*  must  be  positive,  so  that  two  real  roots  exist.  In  formula 
(3.5)  for  9*,  the  quantity  under  the  square  root  h°s  magnitude  less 
than  | 2 oc?  - G|  , since  8c  (r  + be!)  Is  negative.  The  expression 
2c£  - b is  positive  since  c.  > 0,  <?  > 0,  b < 0,  and,  therefore,  hot  . 
roots  are  positive.  From  formula  (3*6)  for  9* , we  note  that  the 
quantity  under  the  square  root  exceeds  | b + 2cd]  In  magnitude. 
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Case  A-i:  B ♦ 2cB  > 0. 

The  value  of  9*  corresponding  to  the  positive  square  root 
satisfies: 


e*  > (aSfr  4 saft).  J 

+ 4C  ’ 

and  hence  is  unacceptable. 

The  value  of  9*  for  the  negative  square  root  satisfies: 

0.  < t sti) . . L 

2c 

A AA  A 

Since  b + 2cd  is  positive,  and  c is  positive,  it  follows 
that  d > -(b/2c),  and  9*  is  an  acceptable  choice  for  a minimum. 

Case  A-ii:  b + 2cd  < 0 (i.e.,  IB  + 2c&l  = -B  - 2cch . 

Here,  the  root  corresponding  to  the  positive  square  root 
satisfies: 


e#  > feed  - b - b - 2c$) 
+ 4c 


b_ 

A 

2c 


Since  B + 2cct  < 0,  d must  be  less  than  -(B/2c);  thus, 
6*  > &,  and  is  unacceptable. 
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Ihe  solution  corresponding  to  the  negative  square  root 


satisfies: 


©*  < 


and  is  acceptable. 


Case  B:  c < 0 (the  coefficient  of  the  quadratic  term  in  fft  is. 

negatived . 

The  quantity  under  the  square  root  must  be  positive  since 

A AA 

c < 0,  r + bd  < 0,  and  thus  there  are  two  real  roots.  In  formula 
(3*5)  for  0*,  the  quantity  under  the  square  root  exceeds  2cct  - b 
in  magnitude,  and  hence  there  must  be  one  positive  and  one  negative 
root.  Because  of  the  sign  reversal  caused  by  division  by  c,  the 
positive  square  root  corresponds  to  a value  6*  < 0,  and  car.  be 
eliminated  from  consideration.  The  only  positive  0*  can  be  written 
as : 


The  quantity  under  the  square  root  has  magnitude  less  than 
| 6 + 2c£|  , so  that  9*  satisfies: 


0*  < 


£ - 2c<?  - b - 2cd 


d. 
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For  Cases  A arvd  B,  then,  assuming  that  f£(Q)  < 0,  the  solution 
e*  corresponding  to  the  negative  square  root  in  formula  (,.5)  or 
(3.6)  will  be  taken  as  the  minimum  of  f^. 

In  order  to  specify  the  function  f^,  four  independent  pieces 
of  information  are  required  about  the  behavior  of  the  function  to  be 
fitted  along  the  search  direction]  f,.  will  therefore  be  used  in 
circumstances  where  a cubic  fit  would  normally  be  carried  out  to 
minimize  a general  function.  Figure  10  illustrates  a barrier  function 
and  its  approximations  by  the  function  f^  and  a cubic,  both  using 
the  same  data.  Hie  more  accurate  modeling  of  the  barrier  function 
by  the  special  function  is  quite  noticeable. 


The  special  function  f^  will  be  fitted  with  the  same  set  of 
data  used  by  a typical  cubic  line  search,  namely,  two  function  values 
and  the  corresponding  two  gradients.  As  in  the  previous  discussion, 
f^  denotes  the  function  value  at  9^,  and  g^  is  the  corresponding 
value  of  the  projected  gradient. 

The  four  unknown  parameters  of  f^  - a,  b,  c,  and  5 — must 
be  solved  for  in  terns  of  the  known  values.  Assuming  that  > 9^, 

the  equations  specifying  f are: 

ft 


There  is  no  loss  of  generality  in  assuming  that  ©^  = 0,  so 
that  from  (A)  and  (B)  tre  obtain  expressions  for  the  coefficients 


a and  $ in  terms  of  d: 


a - f^  + r £n(£) 


b ■ Si  - 5 


Substituting  for  b in  equation  (D) , we  obtain  an  expression 

A * 

for  the  parameter  c in  terns  of  d: 


(e2  - ei + - ~ — ) • 


d d - 0, 


These  three  expressions  for  a,  b,  and  c can  then  be  substituted 
into  equation  (c),  to  obtain  the  following: 


(5-7)  in(cf  - ©J  - in  (3)  + ^ ^ ” — 

d d - 


" r (i^  ^gl  + g2^  ” " fl^  * 


If  a value  d,  the  location  of  the  singularity,  can  be 
found  which  satisfies  (5*7),  the  values  a,  b,  c,  and  9*  can  then 
be  computed  from  the  previously  derived  relationships. 


51' 


wp 


Consider  (3.7) 


as  a nonlinear  equation  in  terms  of  the  function 


where 


(d)  = in 


(^) 


!a(i 

2 \d 


+ g2^  " ^2  ~ 


A 

The  problem  to  be  solved  then  becomes  that  of  finding  a solution  d 
to  satisfy 

(3-8)  *3(d)  = 0. 

If  we  introduce  the  variable  z = 1 - 9g/d  « (d  - e^/d, 
where  0 < z < 1,  we  can  then  write  the  equation  $j(d)  * 0 in 
terms  of  the  variable  z as: 


I 


I 

} 


(3.9) 


in(z)  + | (~  - z)  = . 


The  nonlinear  equation  (3«9)  could  be  solved  for  a suitable 
z,  but  the  function  represented  is  extremely  ill-behaved.  As  z -*  0, 
the  logarithm  term  is  approaching  (-«•),  while  the  reciprocal  term  is 
simultaneously  approaching  (+®)j  it  is  evident  that  (3-9)  is  quite 

fc 

unsuitable  for  purposes  of  computation.  w 

-* 

.1 

■3 
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However,  the  relationship  (3«9)  can  he  transformed  into  an 
equivalent  form  that  is  computationally  reasonable.  A further 
change  of  variable  is  made: 


v = in(z)  , 
v X «v 

so  that  z = e , and  ^ - e $ note  that  v will  be  nonpositive 
since  0 < z < 1. 

When  written  in  terms  of  v,  the  relation  (3.9)  becomes: 
v + \ (e“v  - ev)  = k y 

Since 

sinh(v)  = | (ev  - e~v), 


the  final  result  is: 


v - sinh(v)  = ky 

I 

The  value  v that  satisfies  this  relation  is  a zero  of  the  function 

T , where 
3 


’f.^(v)  * kj  + sinh(v)  - v. 
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Differentiating  with  respect  to  v,  we  obtain: 

» cosh(v)  - 1 
= sinh(v), 

so  that  the  function  has  the  following  properties: 
^j(O)  = ( v = 0 corresponds  to  d »)  j 

lim  *,<v)  * - a; 


> 0 for  v < 0> 

*5(0)  = 0; 

lim  (v)  = »; 
v ~+  - 00  ^ 


<0  for  v < 0. 

These  conditions  imply  that  has  a unique  zero  in  (-  *>,  0) 
if  > 0.  This  requirement  means  that  the  average  of  the  gradients 
at  and  ©2  must  exceed  the  slope  of  the  straight  line  Joining 
f^  and  fg.  If  the  function  to  be  approximated  were  quadratic,  the 
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average  of  the  slopes  at  0^  and  9g  v-ould  exactly  equal  the  slope 
of  the  line  Joining  f^  and  f£.  The  condition  k^  > 0 thus 
implies  that  the  function  to  be  approximated  is  rising  more  rapidly 
than  a quadratic. 

Figure  11  illustrates  the  behavior  of  Vy 

Although  the  function  ¥ unbranded  below  as  v -♦  - » 

(i.e.,  when  z -»  0,  or  d -*  +),  this  property  does  not  cause  any 

computational  difficulties  in  the  current  context.  The  unbounded 
behavior  of  occurs  when  the  estimated  value  of  the  singularity 
is  very  close  to  if  a tolerance,  say  e,  is  specified  such  that 

any  estimate  of  $ is  required  to  satisfy  3 > 9^/ (1  - «) , then  the 
variable  z is  bounded  below  by  €,  and  the  variable  v is  bounded 
below  by  -M,  M > 0,  where  M = - ln(e).  If  values  of  v are 
restricted  to  the  range  (-M,  0) , the  region  where  is  unbounded 

is  eliminated.  If  is  very  large,  it  is  possible  that  the  value 

T^(-M)  will  not  be  negative  for  the  particular  value  of  M chosen, 
and  hence  no  zero  of  will  exist  in  (4/1,  0).  Under  these 

A 

circumstances,  we  simply  accept  v = -M  as  the  solution,  so  that 
£ - 1 - c). 

We  can  easily  solve  the  equation  T^(v)  = 0 with  Newton's 
method,  considering  the  following  properties  of  ¥ . Since  '**'  < 0 
throughout  the  interval  of  interest,  if  the  starting  point  is  chosen 
so  that  < 0,  the  Newton  iterates  will  undershoot  the  solution,  and 
cannot  diverge.  Furthermore,  the  condition  *_(v)  = 0 can  be  written  as: 
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kj  = v - sinh(v) , 


t _ 


and  the  expression .on  the  right-hand  side  is  independent  of  the 
problem  data.  Hence,  the  function  v - sinh(v)  can  be  tabulated 
for  v in  the  range  (-M,  0),  and  by  table  lookup  using  the  value 
ky  a highly  accurate  initial  estimate  of  v,  with  < o,  can  be 
obtained. 

The  Newton  procedure  converges  extremely  rapidly  except  when 
the  value  of  v is  close  to  zero,  because  ¥^(0)  = 0.  This  situation 
is  quite  unlikely  because  the  estimated  singularity  would  then  be 
much  larger  than  For  completeness,  however,  we  note  that  the 

problem  can  be  solved  successfully  even  for  very  small  v.  Hie 
function  ^(v)  can  be  written: 


^(v)  = k^  + sinh(v)  - v 


+ 2 (®V  “ e"V)  " v 


1 , v£  3 

k5  + 2 d + V 4 £-  + r + 


- (1  - V + + •••  ))  -V 


k + — ( 2m  + ^ ■—  + • « • ) - 

+ o + g + 120  ; 


2 


+ + o(v5). 
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For  small  v , the  equation  ¥^(v)  » 0 thus  essentially  becomes 
the  condition 


k 


3 


- 0, 


l/3 

with  explicit  solution  v*  - (-6kj)  ' . Because  we  are  ignoring 
negative  higher-order  terms,  this  value  v*  will  be  to  the  left 
of  the  correct  v,  and  Newton's  method  cannot  diverge.  However, 
the  estimate  v*  is  so  accurate  that  no  iteration  at  all  is  necessary 
to  obtain  an  acceptable  solution. 


4.  iraiaagttkaJUtea 

Hie  safeguarded  linear  searohes  based  on  quadratic  or  cubic 
interpolation  (of.  dill  and  Murrey,  1974)  have  been  modified  for 
use  with  the  logarithmic  barrier  function  by  allowing  interpolation 
with  the  special  functions  described.  Several  rather  complicated 
modifications  are  required  in  order  to  create  an  efficient  algorithm. 
If  no  constraint  is  decreasing  along  the  current  search  direction, 
or  if  no  constraint  approaches  zero  until  sufficiently  far  beyond 
the  starting  point,  then  the  singularity  introduced  to  preserve 
feasibility  will  have  no  significant  effect  on  the  location  of  the 
minimum,  and  the  usyal  linear  search  procedure  should  be  followed. 
There  is  no  computationally  reasonable  way  i,o  determine  & priori 
whether  these  conditions  exist  because  the  constraints  and  objective 
function  may  be  highly  nonlinear,  and  the  effort  expended  to  compare 


hi 
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the  location  of  the  nearest  constraint  zero  with  the  prediction  of 
the  barrier  function's  minimum  might  be  better  used  directly  to 
minimize  the  barrier  function.  The  procedure  to  be  described  seems 
to  be  a satisfactory  compromise  between  excessive  safeguards  and 
unwarranted  assumptions  of  linearity  or  smoothness. 


rhe  choice  of  the  first  step  along  the  search  direction  at 
whioh  the  function  is  to  be  evaluated  is  affected  by  the  possibility 
that  a constraint  may  become  nonpositive  if  the  usual  choice  of  step 
for  the  algorithm  is  taken.  For  example,  with  a Newton- type  method, 
the  initial  step  taken  along  the  searoh  direction  is  unity)  for  a 
quasi-Newton  method,  there  is  normally  a procedure  associated  with 
the  method  for  choosing  the  initial  step.  Let  ccu  denote  the 
initial  step  that  would  be  taken  for  a particular  unconstrained 
algorithm  if  used  to  minimize  the  barrier  function  along  the  given 

A 

direction.  If  a constraint  might  become  zero  at  a < a:  , clearly  a 
shorter  step  than  0!^  should  be  taken.  One  possible  method  for 
determining  the  initial  step  is  to  find  a highly  accurate  estimate 
of  the  step  to  the  nearest  zero  of  a constraint,  say  a,  and  test 
whether  a < Q^.  A subroutine  is  available  that  will,  with  high 
reliability,  locate  the  zero  of  the  nearest  constraint  by  use  of  a 
combination  of  safeguarded  zero-finding  techniques.  However,  locating 

A 

a generally  requires  several  constraint  evaluations,  and  it  may 
turn  out  that  a exceeds  or  is  very  close  to  a , so  that  these 


- a . . u.’i  vUi ifJi ~ib9Jm£bh  i ‘-iUeh  ;£•  ti . •.  uu  uszi 


evaluations  were  essentially  redundant.  One  might  think  of  using 
the  zero-finding  technique  until  the  zero  has  been  shown  conclusively 
to  lie  beyond  a , but  this  approach  involves  quite  complicated  house- 
keeping, and,  more  significantly,  may  still  require  constraint  calcu- 
lations that  do  not  advance  the  computation. 

With  the  "compromise"  algorithm,  the  initial  step  to  be 

taken  along  the  search  direction,  p,  is  computed  as  follows: 


1.  Compute  a , the  step  normally  taken  by  the  unconstrained  method} 

2.  Compute  the  gradient  of  each  constraint  along  p,  i.e., 

T 

8jp,  where  a^  is  the  gradient  of  c^.  For  all  i such  that 
this  gradient  is  negative,  i.e.,  the  i-th  constraint  is  locally 
decreasing  along  p,  compute  = -c^/a^p,  the  predicted  Newton 
step  to  the  zero  of  o^.  Find  5 * min  (c^),  snd  let  I be  the 
index  for  which  a -•  aj.  In  other  words,  3 is  the  smallest 
positive  first-order  step  to  a constraint  zero.  If  no  constraint 
is  decreasing  along  p,  set  the  initial  step  a ^ to  a , and 
skip  the  remaining  logic. 

5,  Estimate  o^,  the  step  to  the  minimum  of  the  barrier  function 
along  p,  which  satisfies: 


(4.1) 


pT(g(x 


V) 


i 

- z 

i-1 


ci(x 


V) 


a^x 


+ V>) 


where  g(*)  is  the  gradient  of  F.  The  relationship  (4.1) 
can  be  used  to  obtain  a crude  estimate  of  ct^,  if  two  approxi- 
mations are  made.  First,  we  assume  that  only  the  Influence 
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....  


of  the  constraint  I is  significant  in  the  location  of 
this  assumption  is  based  on  the  idea  that  for  small  r,  only 
4He  singularity  along  p closest  to  the  starting  point  effects 
the  local  behavior  of  the  barrier  function.  Second,  it  is 
assumed  that  the  gradients  of  the  objective  funotion  and  I-th 
constraint  remain  fixed  locally.  Under  these  conditions  — 
ignoring  second-order  terms,  and  all  but  the  I-th  constraint  — 
(4.1)  becomes: 


c 


T 

r aTP 

i + “t  8ip 


0 , 


so  that 


°b 


°r 

_E_  _ Jli 

1 T 

P ®jP 


Since 


1 


I 


the  estimate  of 


% 


is  given  by: 


«b 


+ 5 . 
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If  p g > 0,  i.e.,  p is  not  a descent  direction  for  the  objective 
function,  then  the  assumption  that  only  one  constraint  will  have 
an  effect,  confined  to  the  neighborhood  of  the  singularity,  is 
not  Justified,  and  we  set  0!^  » y5,  where  0 < y < 1. 

4.  If  a exceeds  Of  , let  a ^ be  as  otherwise,  = a . 

D U U ^ 

In  this  way,  the  initial  step  can  be  chosen  based  on  the 
estimated  decrease  of  the  constraints  if  it  seems  that  their  effect 
will  be  significant  in  the  location  of  the  minimum  of  the  barrier 
function. 

4.2.  jj^MlkUiiaL jSbasJs. 

For  each  step  to  be  taken,  the  set  of  constraints  is  evaluated 
in  order  to  assure  that  feasibility  is  never  violated. 

(k) 

If  any  constraint,  say  the  J-th,  is  non-positive  at  x + cr  'p, 
the  secant  step  to  the  predicted  zero  of  that  constraint  is  computed, 
and  the  next  estimate  of  the  location  of  the  minimum  is  computed  as 
described  above  in  Step  (3),  where  the  secant  step  5j,  given  by: 

- qffi  ef  (x) 

^ c | (x  + aP^p)  - ’ 

is  used  as  a.  This  procedure  is  subject  to  the  safeguarding  require- 
ment that  the  constraints  not  be  evaluated  at  points  that  are  too 
close  together  (see  Gill  and  Murray,  1974)  Brent,  1973,  for  a discus- 
sion of  this  aspect  of  safeguarded  linear  searches). 
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k.3.  Normal  Iteration  of  Modified  Linear  Search 


The  special  functions  are  fitted  during  the  iteration  if  a 
flag  has  been  set  to  'true'.  The  flag  is  set:  (1)  when  the  initial 

was  a^,  implying  that  the  influence  of  some  constraint  is 
predicted  to  be  significant  in  locating  the  minimum)  and  (2)  when 
any  negative  constraint  value  is  encountered  during  execution  of  the 
linear  search,  since  it  has  then  been  demonstrated  that  the  current 
iterate  is  in  a region  influenced  by  the  singularities. 

The  special  functions  are  fitted  iteratively,  .using  the  same 
criteria  for  replacing  points  as  the  usual  linear  searches,  until  the 
particular  convergence  criteria  are  satisfied.  There  are  a few 
subtle  difficulties  in  that  for  small  r,  it  may  be  difficult  to 
locate  the  minimum  because  the  distance  from  the  singularity  to  the 
minimum  may  be  less  than  the  spacing  required  for  constraint  evalua- 
tions. However,  a careful  regulation  of  the  tolerances  involved,  so 
that  impossible  accuracy  is  not  sought,  will  assure  that  the  process 
will  work  as  desired. 

Comparison  with  Usual  Linear  Searches 
In  order  to  determine  whether  the  special  linear  searches  are 
worthwhile,  numerical  experiments  were  oarried  out  for  several 
barrier  functions,  with  varying  values  of  r,  the  barrier  parameter, 
and  t),  the  linear  search  convergence  parameter.  For  the  cubic 
case  where  gradients  are  evaluated  at  every  point,  the  linear  search 
usually  terminates  when 


146 


i .twi'v  Mi  gg-uivu viJi  f 1 hiaiflk t f niii Bhi nmr; uir uftfln 'm*ki  t amv 


I g(x  + at>)|  < n I I(x)| , 


where  g(.)  is  the  projected  gradient  of  the  function  to  be  mini- 
mized) in  the  quadratic  case,  the  linear  search  is  usually  terminated 
when  the  minimum  of  F is  known  to  be  bracketed  in  the  interval 

la,  b J and 


I rta-t.taL-T-Ett-i-aE).  | <,||(l)|, 

i.e.,  when  the  linearised  approximation  to  the  gradient  at  x + oj> 
satisfies  the  same  test  as  g(x  + cqp)  in  the  cubic  case.  There 
are  other  occasions  when  the  normal  linear  search  procedure  will  ter- 
minate, involving  sufficient  smallness  of  the  interval  of  uncertainty, 
closeness  to  the  maximum  permitted  step,  etc. 

Numerous  runs  (about  10)  were  made.  For  both  the  special 
functions  and  the  usual  polynomials,  the  same  initial  step  was  taken, 
and  the  same  procedure  was  followed  for  determining  the  next  point 
if  a constraint  became  negative  during  the  line  search  iteration. 

Hence,  the  only  difference  was  in  the  use  of  the  minima  of  special 
functions,  rather  than  of  cubic  or  quadratic  polynomials,  to  yield 
the  next  point  at  which  the  function  and  constraints  are  to  be  evalu- 
ated as  the  linear  search  proceeds.  In  every  case,  use  of  a special 
function  reduced  the  number  of  function  and  constraint  evaluations; 
the  reduction  became  progressively  more  significant  as  the  value  of 
v was  reduced.  The  reduction  in  the  number  of  function  and  constraint 
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evaluations  ranged  from  7$  to  20$  when  a Merton-type  algorithm  was 
used  to  carry  out  the  unconstrained  minimisation,  arvi  from  12$  to 
2U$  for  a quasi-Newton  algorithm;  thus,  there  was  clear  improvement 
with  the  approximation  by  special  functions. 

5.  tencluaigng 

lhe  extra  work  required  in  the  linear  search  procedure  to  fit 
these  special  functions  is  small.  Some  of  the  housekeeping  (checking 
for  feasibility,  etc.)  must  be  carried  out  with  barrier  functions 
regardless  of  whether  special  functions  are  used  or  not.  The  formu- 
lations presented  here  allow  calculation  of  the  minimum  of  the  fitted 
functions  with  the  same  information  required  to  fit  the  usual  poly- 
nomials. The  singularity  must  be  located  through  an  iteration,  but 
because  of  the  special  form  of  the  iteration  functions,  we  are  able 
to  obtain  a highly  accurate  starting  guess;  in  fact,  in  two  of  the 
three  cases,  the  solution  could  be  obtained  from  tables.  The  itera- 
tion functions  are  well-behaved,  and  Newton's  method  will  usually 
converge  to  the  desired  accuracy  within  two  iterations.  Each  itera- 
tion to  locate  the  singularity  requires  evaluation  of  a transcendental 
function,  but  the  subsequent  reduction  in  the  number  of  function  and 
constraint  evaluations  required  to  locate  a satisfactory  approximation 
to  the  minimum  of  the  barrier  function  seems  ample  justification  for 
use  of  the  special  linear  searches  designed  to  minimize  barrier  functions. 
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Linear  search  algorithms  are  developed  for  use  when  minimiz- 
ing logarithmic  harrier  functions,  whose  one -dimensional  behavior 
is  *n  general  modeled  poorly  by  the  low-order  polynomial  approxima- 
tions of  standard  linear  search  procedures.  Hie  new  methods  are 
based  on  rpecial  approximating  functions  with  a logarithmic  singula- 
rity, and  are  designed  to  utilize  the  same  information  as  procedures 
based  on  special  approximating  functions  with  a logarithmic  singula- 
rity, and  are  designed  to  utilize  the  same  information  as  procedures 
based  on  quadratic  or  cubic  polynomials.  Although  the  parameters 
of  the  fecial  approximating  functions  depend  nonlinearly  on  the 
available  data,  the  determination  of  the  parameters  requires  little 
additional  work  in  comparison  with  polynomial  fits.  Use  of  the 
special  approximating  functions  has  led  to  a significant  improve- 
ment in  efficiency  when  minimizing  logarithmic  barrier  functions, 
where  efficiency  is  measured  by  the  number  of  function  (or  function 
and  gradient)  evaulations  required  for  termination  of  each  linear 
search. 
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